function [F21Ord1, IRF_scaled, structural_shocks]= vm_mult_simple(F,J,Hor,eeta,U,T,n,p,nP,GDPRATIO)
% function [F21Ord1, IRF_scaled, structural_shocks]= vm_mult_simple(F,J,Hor,eeta,U,T,n,p,nP,GDPRATIO)
% Inputs:
%  - F: companion form
%  - J: shock linear combination selector matrix
%  - Hor: impulse response horizon
%  - eeta: elasticities matrix of fiscal with respect to fiscal matrices
%  - U: matrix of VAR residuals
%  - T: number of time periods
%  - n: number of variables
%  - p: number of lags
%  - nP: number of fiscal variables
%  - GDPRATIO: share of non-fiscal variables in GDP
%
% Outputs:
%  - F21Ord1                [Hor by 2 double]   GDP multipliers of G (first column) and T
%                                               shock (second)
%  - IRF_scaled             [Hor by nvars by 3 double]   IRFs to of G shock (first column) and T
%                                               shock (second), and G shock with opposite ordering (third)
%  - structural_shocks      [T by 2 double]     matrix with G (first column) and T 
%                                               shocks (second)

% Based on replication codes for Dario Caldara and Christophe Kamps (2017), The Analytics of SVARs: A
% Unified Framework to Measure Fiscal Multipliers, Review of Economic Studies (2017) 84, 1015–1040
%
% Copyright: 2017 Dario Caldara and Christophe Kamps
% Copyright: 2020-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer

U1 = U(:,1:nP); %fiscal block
U2 = U(:, nP+1:end); %non-fiscal block
u1hat = U1 - U2*eeta; %eliminate feedback (implies partialling out)
b11iSig1 = (u1hat'*u1hat)\u1hat'*U1; %IV regression coefficient of fiscal reduced form residuals on fiscal shocks
b21iSig1 = (u1hat'*u1hat)\u1hat'*U2; %IV regression coefficient of nonfiscal reduced form residuals on fiscal shocks
S1S1 = (u1hat'*u1hat)/(T-n*p-1); %covariance matrix of fiscal shocks
%% deal with correlation between fiscal instruments
% Compute two Cholesky orderings for T and G
s1 = sqrt(S1S1(1,1));
a  = S1S1(2,1)/s1;
s2 = sqrt(S1S1(2,2)-a^2);
SigmaOrd1 = [s1 0 ; a s2]; %T ordered first
s2 = sqrt(S1S1(2,2));
b  = S1S1(1,2)/s2;
s1 = sqrt(S1S1(1,1)-b^2);
SigmaOrd2 = [s1 b; 0 s2]; %G ordered first
ffactorOrd1 = [b11iSig1 b21iSig1]'*SigmaOrd1;
ffactorOrd2 = [b11iSig1 b21iSig1]'*SigmaOrd2;
% first ordering
OmegaTemp  = [ffactorOrd1;zeros((p-1)*n,size(ffactorOrd1,2))];
LtildeOrd1 = vm_irf(F,J,ffactorOrd1,Hor,n,OmegaTemp);
% second ordering
OmegaTemp  = [ffactorOrd2;zeros((p-1)*n,size(ffactorOrd2,2))];
LtildeOrd2 = vm_irf(F,J,ffactorOrd2,Hor,n,OmegaTemp);

%get structural shocks
Q=chol((U'*U)/(T-n*p-1),'lower')\ffactorOrd1; %recover first two columns of rotations matrix from ffactorOrd1=L*Q
structural_shocks=(Q'*(chol(cov(U),'lower')\U'))'; %recover structural shocks from Q'*L^(-1)*u=epsilon where Q'=Q^(-1)

F21Ord1(:,1) = -squeeze(LtildeOrd1(:,3,1))/SigmaOrd1(1,1)/GDPRATIO(2); %get GDP-T multipler, 3rd variable and first shock
F21Ord1(:,2) =  squeeze(LtildeOrd2(:,3,2))/SigmaOrd2(2,2)/GDPRATIO(3); %get GDP-G multipler, 3rd variable and second shock

IRF_scaled(:,:,1) = -squeeze(LtildeOrd1(:,:,1))/SigmaOrd1(1,1)/GDPRATIO(2);
IRF_scaled(:,:,2) = squeeze(LtildeOrd1(:,:,2))/SigmaOrd1(2,2)/GDPRATIO(3);
IRF_scaled(:,:,3) = squeeze(LtildeOrd2(:,:,2))/SigmaOrd2(2,2)/GDPRATIO(3); %get opposite ordering
